Coherent regimes of globally coupled dynamical systems 
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The paper presents a method by which the mean field dynamics of a population of dynamical 
systems with parameter diversity and global coupling can be described in terms of a few macroscopic 
degrees of freedom. The method applies to populations of any size and functional form in the region 
of coherence. It requires linear variation or a narrow distribution for the dispersed parameter. Al- 
though being an approximation, the method allows us to quantitatively study the collective regimes 
that arise as a result of diversity and coupling and to interpret the transitions among these regimes 
as bifurcations of the effective macroscopic degrees of freedom. To illustrate, the phenomenon of 
oscillator death and the route to full locking are examined for chaotic oscillators with time scale 
mismatch. 
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Introduction. Populations of globally coupled dynam- 
ical systems represent a useful framework to study the 
collective properties of biological systems Q, allosteri- 
cally activated enzymatic reactions [2| , electronic devices 
0, and chemical reactions 0, 0, |g7 Although started 
as qualitative analyses, in the last years it has appeared 
that such theoretical approaches can be linked quantita- 
tively to experimental systems, like arrays of Josephson 
junctions |jj or metabolic synchronization in suspensions 
of yeast cells [U • Globally coupled systems can show dif- 
ferent macroscopic behavior when the coupling strength 
and the parameter dispersion are changed. If the cou- 
pling is high enough and the parameter dispersion suf- 
ficiently small, the elements of the population evolve in 
time close to each other (and thus to the mean field) in 
phase space. In the opposite case, the elements of the 
population move incoherently and eventually their posi- 
tions average out, so that the asymptotic dynamics of 
the mean field is characterized only by fluctuations that 
vanish in the limit of an infinite number of elements. Be- 
tween these two limiting cases, complex collective be- 
havior arise. This type of scenario has been observed 
in a wide number of systems, starting from the early 
works of Winfree and Kuramoto on phase rotators IjJ, 
over limit cycle oscillators 01 to chaotic oscillators • 
Investigations have been performed with various statisti- 
cal methods, based on phase reduction 0, the con- 
tinuity equation formalism |13|. or the slaving principle 
(l4|. In particular, the introduction of order parameters 
has appeared to be useful for quantifying the collective 
regimes, as it aims at providing a direct link between 
the microscopic and the macroscopic dynamics. How- 
ever, due to the fact that for low coupling the system 
spans a region of phase space whose dimension grows 
with the population size, the onset of macroscopic oscil- 
lations from incoherence cannot, in general, be described 
by means of a closed system of a few macroscopic de- 
grees of freedom. In the coherent region, on the other 
hand, the macroscopic dynamics lives on a low dimen- 



sional manifold which coincides, in the limit of identical 
oscillators, with the trajectory of a single uncoupled ele- 
ment of the population, and thus a description in terms of 
a few effective degrees of freedom may work. Indeed, for 
limit cycle oscillators with strong coupling and small pa- 
rameter diversity, an order parameter expansion was re- 
cently given by De Monte and d'Ovidio 1 1 5| . showing that 
the transient and asymptotic dynamics of the mean field 
can be accounted for by only two macroscopic variables. 
In the present work we propose a method by which the 
mean field dynamics of globally coupled dynamical sys- 
tems (maps or ODEs) with dispersion on one parameter 
can be systematically reduced (in the coherent regimes) 
to a system of two coupled order parameters. The condi- 
tions we need to impose are either a linear dependence or 
a narrow distribution for the dispersed parameter and, in 
general, smooth dynamical systems. The method can be 
applied to both finite and infinite population sizes. The 
paper is organized in two parts. First, we generalize the 
method proposed in Ref. [15| to any population of ODEs 
or maps: the equations of motion of the mean field and of 
a second order parameter, appearing when the parameter 
diversity is introduced, are derived from the microscopic 
equations through a perturbative approach. The theory 
provides approximate but still quantitative predictions, 
which are tested on populations of ODEs with time scale 
mismatch. In particular, we study the phenomenon of 
oscillator death, in which an attracting equilibrium ap- 
pears in the system due to the interplay between high 
parameter dispersion and strong coupling. We then show 
that the approach continues to apply to non-stationary 
regimes: varying the parameter spread, the transition 
from oscillator death to full synchronization in a popu- 
lation of Rossler systems is followed, and this allows us 
to characterize the various locked regimes in terms of a 
cascade of macroscopic bifurcations. 
Order parameter expansion. Let us consider a popula- 
tion of N dynamical systems in the coherence region, i.e., 
such that the distance \ \X — x-j\\ between the mean field 
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X and the state Xj of any of the oscillators remains small 
in time. For ODEs the dynamics of the j-th element of 
the population is defined by the equation: 

^=f( Xj , Pj )+K(X), 

while for maps: 

Xj -> f(x j ,p j ) + K(X). 

To simplify the notation, we will write in general Txj = 
f( X j,pj) + K(X) where Txj is the operator of time evo- 
lution applied to Xj. The variables Xj G W 1 are the 
state vectors and pj £ 1 is a real parameter taken 
from a distribution of average po = (p) and variance 
<j2 = {{p ~ Po) 2 )- The oscillators are assumed to be glob- 
ally coupled to the mean field X — (x) = Ylj=i X j/N 
through the function K : M™ — > M™. A change of vari- 
ables xj = X + ej expresses the position of every element 
of the population in terms of the mean field, so that in 
the coherent regime a series expansion in ej can be per- 
formed. Setting Sj = pj — po and writing the mixed 
differentials in x and p as T> XiP , we get: 

Txj = f(X,p ) + V x f(X,p )e j + V p f(X,p )5 j 

+V x>p f(X,p )ejSj + o(\e 3 \ 2 ) + C. (1) 

The terms in o(|e.,-| 2 ) are small because of the hypoth- 
esis of coherence. C contains all the Taylor terms that 
are not 0(\ej\ 2 ), i.e., terms with a factor of the type: 
T^x,p m 5 m with to > 1. This term may be neglected if 
the dependence on the parameter is linear or if the pa- 
rameter distribution is narrow. In fact, in the first case 
T> x , Pm 8 m = 0, Vto > 1, while in the second the term 
is 0(6j). We can thus discard these two terms and con- 
struct the macroscopic equations, by averaging Eq. Q 
and using ( e } = (5) =0: 

TX = f(X,p ) + V XtP f(X, Pa )(Se). 

The only first order term that is left defines a new macro- 
scopic variable W := (Se), that we will call the shape 
parameter. It measures the dispersion in both the pa- 
rameter and the phase space. 

In order to close the system, the equation of motion for 
W has to be obtained as well. As for the mean field, this 
can be done by writing: TW — T(5e) — (STx) and 
then using Eq. (JJJ. Discarding again higher order terms 
and making the closure assumption that the term ( ed 2 ) 
is negligible relative to ( 5 2 } (which is true in the coher- 
ent regimes), a closed (approximated) system is obtained 
for the two macroscopic variables: 

f TX = f(X, Po ) + V x . p f(X, Po )W + K(X) 

\ (2) 

{TW = <T 2 V p f(X, Po ) + [V x f(X,p )} W. 



The mean field thus behaves like an uncoupled individual 
element if the oscillators are identical (since in this case 
W = 0). However, parameter diversity may induce new 
regimes. We remark that the method can be applied in 
the same way for other choices of the coupling term. 
Populations with time scale diversity. Let us now apply 
the order parameter reduction to the coherent behavior in 
populations of oscillators with time scale diversity. These 
oscillators may be regarded as a generalization of oscil- 
lators with different natural frequencies (e.g., the Ku- 
ramoto model) and are defined by the equations: 

d ^L=r j g{xj)+lC{X-Xj). (3) 

Tj are strictly positive numbers (which rescale the speed 
along the orbit) and the matrix K, — kT provides an 
isotropic diffusive coupling. 

Oscillator death. As a first application, we will address 
the phenomenon of oscillator death, that was first de- 
scribed in limit cycle oscillators with a large natural fre- 
quency mismatch and strong coupling [Til Hi| . Under 
these conditions, all the elements eventually collapse on 
the origin, which is an unstable focus for the uncoupled 
systems. It was then supposed that the phenomenon is 
more general ^3 an d it was actually described in many 
other systems (like Brusselators and biological sys- 
tems where an equilibrium, unstable in the uncou- 
pled case, becomes attracting and suppresses the peri- 
odic or chaotic oscillations. When the oscillator death 
takes place, the regime is trivially coherent since all the 
elements of the population lie on the same equilibrium 
(which does not need to be unique) and the dependence 
on the dispersed parameter is linear, so that the require- 
ments for applying the order parameter reduction are 
both satisfied. The order parameter expansion will al- 
low us to treat the problem in general and with a simple 
analysis, showing that an unstable equilibrium of g (that, 
without loss of generality, we put at the origin ) can 
always be stabilized by large time scale dispersion and 
high coupling, provided that the equilibrium is a saddle- 
focus whose repelling eigenvalues have an imaginary part 
larger than the real one. In the reduced system Eq. J2J, 
the macroscopic equilibrium corresponding to oscillator 
death is the point X = 0, W — and its stability can 
be studied with simple algebra (the Jacobian matrix is 
composed by four blocks which can be simultaneously 
diagonalized). For each eigenvalue A' = a ; + ibi of the 
unstable equilibrium of g{x), two macroscopic eigenval- 
ues appear: 

a;, 2 =a,-! ± ^/(!) 2 + ^ <4) 

It is clear that, if the equilibrium for the uncoupled sys- 
tem is repelling due to a purely real and positive eigen- 
value, then the corresponding A always has a positive real 
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part and the oscillator death solution cannot be stable. 
Instead, if the instability is due to couples of complex 
conjugates the condition for the oscillator death to be 
stable is that the real parts of the eigenvalues are & U 
negative. 

Let us now consider the case in which only two eigen- 
values of the uncoupled system are complex conjugate, 
while the others are real and negative. This is the only 
case where oscillator death can appear in populations of 
globally coupled two- or three-dimensional systems, like 
coupled Lorenz or Rossler oscillators. The bifurcation 
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FIG. 1: Comparison of the bifurcation boundaries for the os- 
cillator death in the full and the reduced system for different 
populations of dynamical systems: Lorenz (p = 10, b = 8/3) 
r = 28 (•), r = 32 (o), r = 50 (*); Rossler (b = .4, c = 8) 
a = 0.01 (A), a = 0.1 (>), a = 0.4 (<). The continuous lines 
are the bifurcation boundary for the corresponding reduced 
systems. Populations of 100 elements and Gaussian parame- 
ter distributions have been used. 

condition Eq. Q now reduces to a single inequality. Re- 
calling the definition of A, the boundary of the oscillator 
death Re (A) = (where the indexes have been dropped 
due to the fact that there is only one pair of complex 
eigenvalues) can be rewritten in the form: 




Here k* = k/a and 7 = b/ a, corresponding to a rescaling 
of the variables so that the real part of the uncoupled 
system's eigenvalues is unitary. The bifurcation bound- 
ary is thus a surface in the space (k* , a, 7). The fact 
that only three parameters account for the stability of 
oscillator death in any population with the aforemen- 
tioned characteristics follows from the fact that oscilla- 
tor death is essentially a local phenomenon, so that the 
nonlinearities influence at most the critical character of 
the macroscopic Hopf bifurcation. In the plane (k*, a), 
the bifurcation boundaries are lines parameterized by 7 
(Fig. nj. Their asymptotic behavior can be computed 
from Eq. J5J and gives k* = (j 2 — 1) a 2 (upper bound- 
ary, to locking) and k* = 2(1 + a) (lower boundary, to 
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FIG. 2: The bifurcation boundaries of oscillator death for dif- 
ferent types of oscillators collapse on the same, analytically 
determined curve (continuous line) after rescaling the param- 
eters k and a. Symbols as in Fig. Q 



incoherence). Moreover, if the imaginary part of A is 
much larger than the real one (i.e., in the limit 7 — > 00) 
a further rescaling is possible: a* = \J ^f 2 — 1 a. In this 
limit, the bifurcation boundaries can be rescaled to the 
same line: k* — a* 2 and k* = 2. Oscillators as different 
as Rossler and Lorenz systems rescale to the same curve 
obtained analytically from Eq. © (Fig. [5J . 
Non-steady coherent regimes. Let us now speculate about 
what happens outside the region of oscillator death. On 
the lower boundary the coupling is weak relative to the 
parameter spread and the system will bifurcate into inco- 
herence, therefore we expect the expansion to no longer 
be valid. On the upper boundary, by contrast, where the 
closure assumptions are satisfied, the phase transition to 
the locked state arises macroscopically as a Hopf bifur- 
cation. The validity of the introduced approximations, 
moreover, is mantained in the locked regime, where the 
oscillators remain close to the mean field. For this reason, 
not only the transition out of the oscillator death, but any 
other bifurcation up to the fully synchronized state can be 
identified through the reduced system. As an example, 
Fig. |3 compares the behavior of the mean field of a pop- 
ulation of 32 Rossler oscillators with time scale mismatch 
to the order parameter expansion Eq. Changing the 
spread a in the time scale distribution, a complete pe- 
riod doubling cascade can be followed, connecting the 
fully locked regime to oscillator death, with remarkable 
quantitative agreement. 

From the validity of the order parameter expansion in 
describing the macroscopic chaotic dynamics one can in- 
fer that the Lyapunov spectrum of the reduced system is 
composed by the dominant eigenvalues of the full system. 
The microscopic degrees of freedom which are neglected 
in the approximation would therefore act as a perturba- 
tion on the collective dynamics and give a small contribu- 
tion to the Kaplan- Yorke dimension of the macroscopic 
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attractor. 

Conclusions. In ensembles of globally coupled oscilla- 
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FIG. 3: Poincare section of the mean field and of the shape 
parameter for a population of 32 coupled Rossler oscillators 
(a = 0.25, b — 1, c — 8.5) with time scale mismatch (top) and 
for its order parameter reduction (bottom). The coupling is 
k = 1. The systems go from chaos to oscillator death when the 
standard deviation a of the parameter distribution increases. 
The reduced system reproduces the bifurcation cascade of the 
population with remarkable quantitative agreement. 

tors, parameter diversity may induce nontrivial collective 
behavior, where the mean field dynamics is qualitatively 
different from that of each uncoupled element. In this 
work we have shown that such regimes are low dimen- 
sional in the region of coherence, and that a descrip- 
tion by means of few effective degrees of freedom may 
be given. This is done through an expansion around the 
perfectly synchronized state. Although the method in- 
volves some approximations (it is not exact even for a 
linear analysis), it provides nevertheless an accurate and 
quantitative description of the dynamics at the macro- 
scopic level. We remark that there are no requirements 
on the population size: as far as the oscillators are of the 
same type, populations of different sizes behave the same 
if they have the same coupling term and the same vari- 
ance of the parameter distribution. The somewhat sur- 
prising consequence of this fact is that the macroscopic 
features of the coherent regimes can be accounted for by 
a system of just two coupled oscillators, providing a two- 
body approximation of the population dynamics. Finite 
size effects, however, do arise as the region of incoher- 
ence is approached, and the dimension of the collective 



dynamics increases. There are several ways in which we 
think our approach can be developed further. In partic- 
ular, the inclusion of higher order terms may allow us to 
explain more complex collective regimes arising close to 
incoherence. Moreover, a similar approach accounts for 
the effect of noise on the collective dynamics of identical 
oscillators. A macroscopic bifurcation scenario similar 
to that induced by parameter diversity appears in large 
populations of noisy chaotic maps. This bifurcations and 
the finite size effects can again be described in the frame- 
work of an order parameter expansion |19| . 
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